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SUMMARY 


The multidimensional , ensemble-averaged, compressible, time-dependent 
Navier-Stokes equations are solved to predict the turbulent flow field 
resulting from confined swirling and nonswirling jets discharging into a 
suddenly expanded duct. The calculations which correspond to the experiments 
of Johnson and Bennett were conducted in a domain whose inflow boundary was 
situated upstream of the dump plane where the flow is unaffected by viscous 
interactions and extended downstream into the duct where the flow is fully 
developed. In order to be faithful to the actual experimental configuration, 
all sharp corners were retained and the inner jet wall was tapered. 

A two-equation k-e turbulence model was employed to obtain the reported 
results. For the swirling case there was excellent qualitative and 
quantitative agreement with the experiments while for the nonswirling case 
qualitative agreement was obtained. The differences in agreement between the 
numerical predictions and the experimental data in the two cases appear to 
correspond to the effect of large scale coherent structures in the flow 
field. As determined in the companion paper by Brondum and Bennett, these 
large scale structures are dominant in the nonswirling case and may thus 
have a significant effect on the turbulence model, thereby leading to the 
discrepancies noted in the numerical computations. Furthermore, the 
calculations show that excessive artificial dissipation can have a dramatic 
effect on the overall flow structure, and must be effectively controlled to 
obtain accurate predictions. 
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INTRODUCTION 


The design of modern gas turbine combusters is an extremely complex 
process. Many factors can influence their operation, including geometric 
effects, inflow properties of the air and fuel, turbulence of the flow and 
the overall mixing process. Obviously, methods that could aid in 
understanding and simulating these phenomena would be extremely useful in 
the design of more efficient combustors. Before considering the combustor 
as a whole includeing the combustion process, it is advantageous to consider 
the isothermal case. This allows study of individual nonreacting fluid 
dynamics phenomena, such as the effects of mass and momentum transport 
prior to considering the complications associated with reacting fluids. 

An understanding of these processes is a necessary prerequisite to the 
simulation of the more complex reacting flow field. 

Recently Johnson and Bennett (Ref. 1) and Roback and Johnson (Ref. 2) 
have accumulated extensive experimental data for nonswirling and swirling 
turbulent flows in confined suddenly expanded coaxial jets. Since the 
experimental configurations were constructed to be similar to actual 
combustors and the Reynolds number under which the experiments were conducted 
was sufficiently high to assure fully turbulent flow, the data obtained from 
these experiments could be directly applied to the study of gas turbine 
combustors. 

References 1 and 2 describe tne non-instrusive experiments conducted on 
a confined suddenly expanded coaxial jet. A schematic of the facililty is 
shown in Fig. 1. The working fluid was water which was circulated by a pump 
from the storage tank through the test section. Laser Velocimeter (LV) and 
Laser Induced Fluorescence (LIF) techniques were employed to obtain the data: 
velocities, concentrations and flow visualization. Details of the operation 
of the system are given in Refs. 1 and 2. 

In association with these experiments numerical computations were 
performed (cf. Ref. 3) to model the flow fields. The results of these 
calculations indicated that several areas were in need of special attention; 
turbulence modeling, specification of upstream (inflow) boundary conditions 
and control of numerical or aritifical dissipation. With regard to 
turbulence modeling it should be noted that, a major objective of the 
experimental program was to obtain data bases from which a better 
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understanding and formulation of transport models could be obtained. 

In accordance with this goal the present effort was initiated to obtain 
accurate numerical computations which is a prerequisite in meeting the aims 
of the program. Concurrently Brondum and Bennett in the companion effort 
(Ref. 4) conducted experiments to isolate the effects of the large scale 
structures. 

The effects of boundary conditions in general and inflow boundary 
conditions in particular play an important role in the flow development and, 
therefore, must be carefully chosen. In view of the strong interactions that 
occur between the coaxial jets, the subsequent strong mixing and the large 
recirculation zones that develop it is generally agreed that the full 
ensemble-averaged Navier-Stokes equations should be considered rather than 
simplified systems of equations. However, there is not universal agreement 
on how and where to specify boundary conditions. Implicit in the choice and 
application of "correct" boundary conditions is the choice of the appropriate 
computational domain. Recent results employing the TEACH code and its 
derivatives consider a rectangular domain with the upstream boundary being 
situated at the dump plane (cf. Ref. 3). In this computational domain one 
must specify at the dump plane not only the streamwise velocity and swirl 
velocity (if swirl Is present), but also the normal velocity which is an 
extremely sensitive quantity. Furthermore, at the dump plane strong 
interactions between the jets occur so that it is not an optimum location at 
which to impose boundary conditions. Since the determiation of the flow 
properties at or near the dump plane is one of the objectives of the 
calculation, one cannot specify boundary conditions there. 

In contrast to the procedure described above in the present calculations 
the upstream inflow boundary is placed upstream of the dump plane, where the 
flow properties are not influenced by the Interaction process occurring at 
the dump plane. Further, the geometrical domain has been constructed to 
include a tapered inner wall for the central jet in order to model as best 
possible the actual experimental facility. 

Navier-Stokes calculations were performed for the two cases considered 
in Refs. 1 and 2 employing a k-e turbulence model. The Navier-Stokes 
solution procedure which was used in this effort was the consistently split 
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linearized block, implicit (LBI) scheme of Briley and McDonald (Refs. 5 and 
and 6). The numerical scheme is embodied in a general computer code termed 
MINT (Multidimensional, Implicit, Nonlinear, Time-Dependent). The particular 
form of the code being used for the present application solves the general 
tensor form of the Navier-Stokes equations and, therefore, can be used with a 
general coordinate system. The dependent variables in the analysis are the 
velocity components, the density, and for turbulent flow if a two-equation 
model is used, the turbulence kinetic energy, k and the dissipation rate, e. 
The results obtained compared well with the experiments in Refs. 1 and 2. 

In particular, for the nonswirling case there was excellent quantitative 
agreement with the data. Furthermore, the calculations indicate that 
numerical dissipation can have a significant effect on the numerical results. 
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ANALYSIS 


The present analysis is based upon the solution of the ensemble-averaged 
Navier-Stokes equations using the linearized block implicit (LBI) method of 
Briley and McDonald (Ref. 5). The equations are solved in a constructive 
coordinate system with density and the velocity components being taken as 
dependent variables. The discussion of the coordinate system and governing 
equations is given next. 


Coordinate System 

In Fig. 2 is shown a schematic of the experimental facility with the 
pertinent dimensions (duplicated from Ref. 2). It consists of two coaxial 
pipes of Inner and outer radii Rij and Ra respectively expanding into a 
larger radius pipe of radius Rq with the inner pipe having a taper angle of 
7.5 degrees. The computational domain chosen for the numerical computation 
was constructed to model this configuration as close as possible and is shown 
in Fig. 3. 

There are several important features that are noteworthy. First, and 
foremost, the computational domain extends upstream of the dump phase, away 
from where the two jets mix. This allows for the specification of boundary 
conditions in a region unaffected by the mixing of the two jets which in 
general cannot be prescribed accurately. However, the computational 
domain introduces additional complications which must be taken into account 
by the solution procedure. These include reentrant corner points at the 
intersection of the Inner pipe walls with the dump plane walls, and the 
introduction of surfaces where boundary layers develop and which must be 
resolved when no-slip boundary conditions are employed. In the case under 
consideration three reentrant corners are introduced. Note that although the 
inner central jet wall is tapered it terminates at the dump plane such that 
there is a finite thickness of the wall between the two jets. Furthermore, 
three additional walls have been introduced at Ri j , Ri 2 and Ra, upstream of 
the dump plane where the respective boundary layers must be resolved. This 
places additional demands on the grid, viz. more grid points and/or grid 
clusterings are required. 

The next feature of note is the inclusion of tapered inner wall at a 
nominal angle of 7.5° which corresponds to the actual experimental setup. 
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This tapered section is not a straight line, but rather consists of a cosine 
curve in order that the coordinate lines vary smoothly. By including this 
tapered portion, the grid becomes nonorthogonal , and must be taken into 
account in the analysis. 

This geometry although more complicated than those considered by other 
researchers (e.g. Ref. 3) is more realistic and is, therefore, employed. 
Nevertheless, the nonorthogonal geometry and the reentrant corners requires 
no additional code modifications since the MINT procedure has been designed 
to treat such cases . 

The upstream boundary was placed at 51.0 mm upstream of the dump plane 
to correspond to the location where the swirler was placed and was employed 
for the nonswirl case as well. The downstream boundary was placed at 14.0 
R 0 downstream, where fully turbulent conditions should be recovered. 

The grid consisted of 91 grid points in the radial direction and 71 
points in the streamwise direction. The grid was clustered as noted above to 
resolve the boundary layers on all solid surfaces. The solid surfaces 
correspond as the following radial grid point locations 

R n = 28 
R i2 = 31 
Ra = 63 
R 0 = 91 

In the streamwise direction, there were 71 grid points with the dump plane 
located at grid point 15. The maximum grid spacing was at the downstream 
boundary, AZ max = .70 Rq, and the minimum spacing at the dump plane was 
AZ m in = .01 Rq. The family of radial lines were straight while the 
streamwise curves conformed to the body shape. As noted above, this led to a 
nonorthogonal coordinate system, and is shown in Fig. 4. Note that 
downstream of the dump plane the streamwise coordiante curves are adjusted to 
cluster more points near the centerline. To achieve the grid clustering, a 
transformation due to Oh (Ref. 8) was employed. 
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Governing Equations 


The equations used in the present effort are the ensemble-averaged, 
time-dependent Navier-Stokes equations which can be written in vector form as 


Cont inuity 


dp — 

d\ r 


( 1 ) 


Momentum 

dpU —— _ 

+ V(pUU) = — VP + V* (ir + tt T ) 

dt 

Energy 


( 2 ) 


dp h — — — D P 

+ V-(pUh) = -V-(Q + Q T ) + +<P + p€ (3) 

dt r Dt r 


where p is density, U ls velocity, p is pressure, tt is the molecular stress 
tensor ir T is the turbulent stress tensor, h is enthalpy, Q is the mean heat 
flux vector, Q T is the turbulent heat flux vector, $ is the mean flow 
dissipation rate and e is the turbulence energy dissipation rate. If the 
flow is assumed as a constant total temperature, the energy equation is 
replaced by 


T = T + = constant 

f 2C n 


(4) 


where T t is the stagnation temperature, q is the magnitude of the velocity 
and Cp is the specific heat at constant pressure. In the cases considered 
in this work, constant total temperature has been assumed. A number 
of terms appearing in Eqs . 1-3 require definition. The stress tensor 
appearing in Eq . 2 is defined as 

TT = 2/xtD -(Y/l- K B )VTJl (5) 

where Kg is the bulk viscosity coefficient, I is the identity tensor, and 
tD is the deformation tensor, defined by: 

id =^-((Vu) + (vIT) T ) < 6 > 
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In addition, the turbulent stress tensor has been modeled using an isotropic 
eddy viscosity such that: 


7 r T = — p u' u‘ = 2fJL T \D - (/t T V- U + pK) I ^ 

where, k is the turbulent kinetic energy and p-p, the turbulent viscosity, 
is determined by a suitable turbulence model. Turbulence modelling is 
described in some detail in the next section. 

Equation 8 contains a mean heat flux vector defined as follows: 

O' = — /c V T <8) 

and a turbulent heat flux vector defined as: 

0 T =-/c T VT (9) 

where < and k t are the mean and turbulent thermal conductivities, 
respectively. 

Also appearing in Eq. 3 is the mean flow dissipation term 

= 2/AtD:0-(^/x-K B )(V-LT) 2 (10) 

The equation of state for a perfect gas 

P = pRT (11) 

where R is the gas constant, the caloric equation of state 

e = C V T (12) 

and the definition of static enthalpy 

h = C p T (13) 

supplement the equations of motion. 


Finally the flow properties p, < and Kg are determined using the following 
constitutive relations. 

The molecular viscosity u is determined using Sutherland's law. 


_ (JL \ 3/2 t q+ s » 

Po = ' V T +S, 


(14) 


where = 10 0°K for air. 
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The balk viscosity is assumed to be zero 


k b = o 


( 15 ) 


and the thermal conductivity is determined by use of a relation similar to 
Sutherland's law viz. 


/ T f/2 T Q + S„ 
if. = WJ T + 


(16) 


where S 2 = 194°k for air. 


Dependent Variables and Coordinate Transformation 


The set of governing partial differential equations which model the 
physical processes was presented in the previous section. For generality 
these equations were written in vector notation; however, before these 
equations can be incorporated into a computer code, a coordinate system must 
be chosen. The governing equations can then be cast in a form reflecting the 
choice of the coordinate system. Therefore, the governing equations written 
in a cylindrical polar coordinate system are transformed with a general 
Jacobian transformation of the form 

y j = y j (x,, x 2 , x 3 ,t) 


T = t 

where (x_, X 2 , X 3 ) are the original coordinates (Ref. 9). In cylindrical 
polar coordinates (x x ) would correspond to (r, 0, z) . The velocity 
components remain the components, ( U 1 , U 2 , U 3 ) in the (xj, X 2 > X 3 ) coordinate 
directions, respectively. The new independent variables yJ are the 
computational coordinates in the transformed system. The coordinate system 
requirements for the problem under consideration may be represented by a 
subset of the general transformation, Eq . (17) 


y' = y l (x l ,x 3 ,n 
y z = y 2 (x 2 ) 
y 3 = y 3 (x,, x 3 , t) 


(18) 
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which is a general axisymmetric time-dependent transformation. For the 
coaxial jet configuration which is axisytranetr ic , Eq. (18) reduces to 
y = X 2 and all derivatives 3/3y are assumed to be zero. 

Application of the Jacobian transformation requires expansion of the 
temporal and spatial derivatives using the chain rule, i.e., 


d4> d (fi 


at 


and 


dr 


3 i 

+ I y,t 

j=l 


y J 

dx ■ j=i ay J 


dyj 


(19) 


( 20 ) 


where 


j _ <3y j 

y,» - 


at 


^ ■ f 


( 21 ) 


The relations Eqs . (19-21) are first substituted into the governing equations 
(1-4) written in Cartesian or cylindrical polar coordinates. Then the 
resulting equations are multiplied by the Jacobian determinant of the inverse 
t rans format ion , 



dX| 

dX| 

ax, 


ay 1 

dy 2 

ay 3 

a(x,,x 2l x 3 ) 

dx 2 

dx z 

dx z 

<3(y' ,y 2 ,y 3 ) 

dy' 

dy 2 

dy 3 


dx 3 

ax ? 

dx 3 


a y ' 

a y 2 

a y 3 


and the equations are cast into a "semi-strong" conservation form (Ref. 9) 
using the following relations, 


3 

I 

j =l 



(23) 
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and 


!£ + y aj y j ,t 

dr j4i ay i 


= 0 


(24) 


The semi-strong conservation form implies that all factors involving the 
radial coordinate r = remain as they were before the Jacobian 
transformation. The resulting equations are presented in Appendix B. 

The geometric relations Eq. (23-24) may be obtained from the 
transformation relations for Jy,j^ and Jy,-} in terms of the inverse 
transformation derivatives (e.g., Ref. 10), 


and 


Jy ,i 


Jy 2 = x 

I 

Jy 2 

Jy 


,2 

,2 

,3 


V t ^ — V 


J y 


Jv 3 


X 2,2 X 3,3 

“ X 2,3 X 3,2 

X 3,2 X i,3 

“ X 3,3 X /,2 

x "x 

1,2 2,3 

- x "x 

1,3 2,2 

X 2,3 X 3,1 

- *2,. *3,3 


‘ X 3,1 *1.3 

V 

- "... "3.3 


'VS,, 

X 3,1 *1,2 


X l,l X 2,2 

~ X 1,2 X 2,l 




y,t 


f 


(25) 


(26) 


11 



Turbulence Modeling 


Several alternative turbulence models can be applied to the problem at 
hand. In general terms, these models are the zero-, and two-equation 
models. The formulation of each of the two is described in this section. 


Zero Equation Model - (Mixing Length) 


Of all available turbulence models, Prandtl's mixing length model is 
probably still the most widely used. The model was originally developed for 
use in unseparated boundary layer flow situations and has been shown to 
perform well under such conditions. An advantage of the method from the 
point of view of economy is that it does not require additional transport 
equations to model the effect of turbulence, but rather relates the Reynolds' 
shear stress to mean flow quantities via: 


I I 




auj 

dXj 


where 

jJLj = />£ 2 (21D : (D) l/2 


and 

£ = min [ /c<Jd] 

where d is the normal distance to the nearest wall and D is the van Driest 
damping coefficient given by 

D = I - exp(-y + /A*) y + = du T /u 

= 0.098 ( 27 ) 

U T - {T l / P ) " Z 

K = 0 4 

and where the local shear stress T£ is obtained from 

r t = <2tD: E>) ,/2 ( 28) 

and OD is defined by Eq. 6. 
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One problem in the mixing Length formulation is the definition of 6, 
which for flows of the type considered here is extremely difficult to 
estimate. Hence, it has not been used to obtain the reported flow fields. 
However, as will be discussed in the following sections, this model is used 
to initialize the computations. 


Two-Equation Model - (k-e) 


As discussed above, the mixing length concept is valid for a variety of 
flows, in which the viscous layer is wall bounded. However, in cases such as 
considered here which involve large reciruclation zones, and is shear 
dominated, a less restrictive model is required. One such model is the 
two-equation turbulence model (Refs. 11-15) in which a transport equation for 
turbulence kinetic energy, k, is formulated as follows: 


dpK 

~dT 


H-r 


1 / 2,2 


+ V-(/3 UK) = V- ( — -VK) + 2/J. t (ID: CD) - p€~2pv{VK u *) 


(29) 


where k is the turbulence kinetic energy and is defined as 

K .-J u*V (30) 

and the transport equation for the dissipation e is 

dp€ + V • (/>U e)= V (^lVe) + C,(2 /x t ID : ID)_L + 2/i/x t ( V 2 U) 2 - C z p _l! (31) 

<3t cr € K K 

However, attempts to solve Eqs . 29 and 31 without modification present 
problems because an appropriate boundary condition for e at a solid boundary 
ls difficult to prescribe such that Eq . 31 is satisfied. Following the 
suggestion of Jones and Launder (Ref. 12), the turbulence dissipation 
equation has been modified by the inclusion of the term: 


-2/i/z T (V 2 U> 2 
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in the energy dissipation equation, Eq . (31), and by the inclusion of the 
term: 

-2 />z/(VK 1/2 ) 2 


in the turbulence energy equation. These additional terms allow an e = o 
wall boundary condition to be applied and appear to correctly model the near 
wall region as discussed in Ref. 12. Following Ref. 12, the following 
empirical relations are used. 

a € - 1.3 cr^ = 1.0 

C } = 1.43 

= 0.09 exp[-2.5/(l + R T /50)j 
C 2 = 1.92 [l.O - 0.3 exp (-R|)] 

and R t is defined as: 


R 


T 


pk 2 


The Prandt 1-Kolmogorov relation, defines the turbulent viscosity as: 


-c P * Z 

H- r- c p.~r- 


(32) 


In modeling the flow in the near wall region where low local turbulence 
Reynolds' numbers occur, two approaches are available. The first is the wall 
function approach which does not resolve the near wall region but assumes 
specific function forms for the required turbulence quantities and uses these 
forms to create the required normal derivative formulations at the first grid 
point from the wall. Such an approach obviously requires a detailed 
knowledge of the turbulence model dependent variables in the vicinity of the 
wall. Although reasonable function formulations can be specified for simple 
two-dimensional flows such as constant pressure boundary layers, 
specification in the much more complex flows of current interest is much more 
difficult. Therefore, the alternative approach in which the viscous sublayer 
is resolved has been used. The method makes no approximation at the 
boundary, but requires that the near wall low turbulence Reynolds' number 
physics be modeled. 
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Initial and Boundary Conditions 


Steady solution of the system of governing partial differential 
equations represented by Equations (1-3) is obtained by time marching these 
equations until a steady state is reached. Before the solution procedure is 
described two important aspects must be discussed: (1) the initial 

conditions and (2) the boundary conditions. Any procedure which utilizes 
either a time marching method to obtain a steady state (or transient) 
solution or a Newton-Raphson iteration procedure requires some initial guess 
of the flow variables (in this case all the dependent variables and other 
necessary variables such as pressure, temperature, viscosity, etc.). In some 
of the simpler cases, some reasonable approximation to a converged solution 
can either be guessed or obtained through physical reasoning. However, since 
the flow field considered under this effort is dominated by large 
recirculation zones it was felt that an initial guess containing such closed 
vortical patterns would be very difficult at best and at worst could hamper 
the ultimate convergence history. The approach taken here was to assume that 
the flow was initially stagnant (all velocity components were set to zero), 
and that the pressure and temperature were constant being set equal to the 
downstream exit flow conditions. The upstream velocity profiles were then 
raised to some precribed level over a period of time thereby driving the flow 
through the duct. Thereafter the solution was marched out in time until a 
steady state was achieved. This technique has the advantage of being easy to 
implement in any geometric configuration. 

To obtain a solution of the governing system of partial differential 
equations represented by Equations (1-3), it is necessary to define boundary 
conditions on each bounding surface of the computational domain. For the 
purposes of this investigation boundary conditions can be classified as 
occurring on two different types of bounding surfaces: (1) walls on solid 

surfaces, (2) inlets and exits. The boundary condition utilized on each 
different type of surface will now be discussed in turn. 

At walls and solid boundaries no slip is prescribed, i.e. the streamwlse 
and normal velocities are set to zero. In addition the normal pressure 
gradient is set to zero. As an alternate boundary condition the normal 
momentum equation can be solved at the boundary, and is employed if required 
by the physics of the flow. 
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At the outflow boundary, for subsonic flow the static pressure is 
specified and the velocity components are extrapolated, i.e. the second 
derivatives of the streamwise and normal velocities are set to zero. The 
inflow boundary, however requires some additional care. Since the velocity 
profiles at the inlet were measured (cf. Ref. 7), they were specified there 
and were fixed throughout the calculation. This essentially sets the mass 
flow through the system. An additional boundary condition is needed for 
density, which for this case reduced to an extrapolation condition. 

As a final note, a description is given of the treatment of the 
reentrant corners, which appear as geometrical singularities in the flow. 

The specification of no slip offers no difficulty at these corners, since the 
velocities are set identically to zero. However, the pressure condition is 
somewhat more difficult since it involves a derivative. In order to 
circumvent the difficulty associated with choosing the direction in which the 
normal derivative is to be taken, the corner is treated as a double valued 
point. Hence, two values of pressure (and density) are stored, each 
corresponding to the direction in which the coordinate lines approach the 
corner. Although this method is approximate, it has worked well in practice 
and does not appear to adversely affect the results obtained. 

Numerical Procedure 

The numerical procedure used to solve the governing equations is the 
consistently split linearized block implicit (LBI) scheme originally 
developed by Briley and McDonald (Ref. 5). A conceptually similar scheme has 
been developed for two-dimensional MHD problems by Lindemuth and Killeen 
(Ref. 16). The procedure is discussed in detail in Refs. 5 and 6. The 
method can be briefly outlined as follows: the governing equations are 

replaced by an implicit time difference approximation, optionally a backward 
difference or Crank-Nicolson scheme. Terms involving nonlinearities at the 
implicit time level are linearized by Taylor expansion in time about the 
solution at the known time level, and spatial difference approximations are 
introduced. The result is a system of multidimensional coupled (but linear) 
difference equations for the dependent variables at the unknown or implicit 
time level. To solve these difference equations, the Douglas-Gunn (Ref. 17) 
procedure for generating alternating-direction implicit (ADI) schemes as 
perturbations of fundamental implicit difference schemes is introduced in its 
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natural extension to systems of partial differential equations. 

This technique leads to systems of coupled linear difference equations having 
narrow block-banded matrix structures which can be solved efficiently by 
standard block-elimination methods. 

The method centers around the use of a formal linearization technique 
adapted for the Integration of initial-value problems. The linearization 
technique, which requires an implicit solution procedure, permits the 
solution of coupled nonlinear equations in one space dimension (to the 
requisite degree of accuracy) by a one-step noniterative scheme. Since no 
iteration is required to compute the solution for a single time step, and 
since only moderate effort is required for solution of the implicit 
difference equations, the method is computationally efficient; this 
efficiency is retained for multidimensional problems by using what might be 
termed block ADI techniques. The method is also economical in terms of 
computer storage, in its present form requiring only two time-levels of 
storage for each dependent variable. Furthermore, the block ADI technique 
reduces multi-dimensional problems to sequences of calculations which are 
one-dimensional in the sense that easily-solved narrow block-banded matrices 
associated with one-dimensional rows of grid points are produced. A more 
detailed discussion of the solution procedure as discussed by Briley, Buggeln 
and McDonald (Ref. 18) is given in the Appendix A. 

Artificial Dissipation 

Since the calculations of interest are often at high Reynolds numbers 
typical of combuster flow fields, it is necessary to suppress spatial 
oscillations associated with central spatial differences approximations. 

This can be done via a dissipative spatial difference formulation 
(e.g., one-sided difference approximations for first derivatives) or by 
explicitly adding an additional dissipative type terra. For the Navier-Stokes 
equations, the present authors favor the latter approach since when an 
additional term is explicitly added, the physical approximation being made is 
clearer than when dissipative mechanisms are contained within numerical 
truncation errors, and further, explicit addition of an artificial 
dissipation terra allows greater control over the amount of non-physical 
dissipation being added. Obviously, the most desirable technique would add 
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only enough dissipative mechanism to suppress oscillations without 
deteriorating solution accuracy. Various methods of adding artificial 
dissipation were investigated in Ref. 20, and these were evaluated in the 
context of a model one-dimensional problem containing a shock with a known 
analytic solution (one-dimensional flow with heat transfer). The methods 
which were considered included second-order dissipation, fourth-order 
dissipation and pressure dissipation techniques. 

As a result of this investigation, it was concluded that a second-order 
anisotropic artificial dissipation formulation suppressed spatial 
oscillations without impacting adversely on acccuracy. In the present 
application a term of the form 




] 


is added to the governing equations for each coordinate direction j . 

The variable <)> denotes the velocity component IJ^ for the -direction 
momentum equation, the density p for the continuity equation, and the 
enthalpy h for the energy equation. The coefficient (p ar t)j is obtained 
from 


/>UjAxj < (i/a d )(/I +(fJi ar t )j) 


where Axj is the grid spacing at the point in question. The quantity p 
denotes the effective viscosity (p e ff) for the momentum equations, 

(Peff/Pr) for the energy equation, (Peff/°k) for the turbulence kinetic 
energy equation, (Peff/ CT e) for the turbulence dissipation equation, and 
is zero for the continuity equation. 

The question arises as to the values of o x and Oy which should be 
chosen. This was assessed both through model problems (Ref. 19), and through 
actual calculations (Refs. 19, 20 and 21). These results indicated that 
values of a - .5 which corresponds to a cell Reynolds number 2 limitation 
would severely damp physical variations. However, when a was set In the 
range .05 < a < 0.10, which correspond to a cell Reynolds number range 
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between 20 and 10, spurious spatial oscillations were damped with no 
significant change In the calculated results as o was varied In this range. 
Further, as discussed in Refs. 19-21, the results obtained showed good 
agreement with data. 
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DISCUSSION OF RESULTS 


As discussed in previous sections, the computation of the flow field was 
initiated from quiescent conditions. During the initial transient phase the 
upstream velocity profile was augmented until its magnitude matched the 
experimental values of Johnson (Ref. 7). Once the profiles were attained, 
they were fixed for the duration of the calculation. 

The velocity profiles that were employed were obtained from curve fits 
of Johnson's experimental data (Ref. 7) which was taken at 41mm upstream of 
the dump plane using a hot film probe. For the inner pipe a log law profile 
in conjunction with a viscous sublayer was fit to the data. With the 
constants in the log law formula fixed, the sole parameter is the friction 
velocity u T , which was determined by a Newton iteration procedure. 

In the annulus, the formulas presented by Bird, Stewart and 
Lightfoot (Ref. 22, Eqs. 5.F-2 and 5.F-3) were used to fit the streamwise 
velocity data. 


u max u 


K, 


r k 2 - X 2 i 1 7 2 r ( X - k ) R I 
U rl k J ln l r-kR I 


< XR 


u max ~ u = X 2 In jy -j r > XR 


where = .4, k = R£/R 0 and r = XRq is the point of maximum velocity. 

Since each formula is singular at the inner and outer walls, these formulas 
were supplemented by additional equations that account for the viscous 
sublayer. Here again the only free parameters were the two friction 
velocities at the respective wails, and were determined by using a Newton 
iteration procedure. 

The k and e distributions at the upstream boundaries were also obtained 
from the experimental data given in Ref. 6. In contrast to the curve fit 
procedure used for the streamwise velocity distribution, for k and C local 

I 

piecewise parabolic polynomials were employed. The curve fits also 
prescribed zero values for k and t at the walls consistent with the boundary 
conditions used in the calculation procedure. 

For the swirling case, the swirl velocity distribution was also required 
at the upstream boundary. Since that data was not given, an alternate 
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procedure was employed to obtain the swirl velocity profiles. A comparison 
of the nonswirling axial velocity profiles at the upstream station and at 5mm 
downstream of the dump plane showed that the core velocity profile away from 
the walls does not vary significantly for the two cases. Hence, it was 
assumed that similar behavior should hold for the swirl or azimuthal velocity 
profile. Therefore, the data given at z = 5mm downstream of the dump plane 
was used to obtain the required velocity profiles. The curve fit procedure 
employed was identical to the nonswirling case in which the data was fit with 
logarithmic profiles in conjunction with a laminar sublayer. 

Both the swirling and nonswirling cases were computed for a Reynolds 
number of 35,000 based on duct diameter. The nonswirling case was run first 
employing an artificial dissipation parameter of a = .5. The results were 
very similar to coarse grid calculations considered previously. Initially, 
the nonswlrl case was computed on a coarser grid of 61 x 51 grid points. 

This calculation differed from that which is described in this report in two 
respects. First, the inner jet wall was untapered and second, the upstream 
velocity profiles were "guessed" since at that time no experimental data was 
available. The results were substantially the same as the current 
nonswirling case which is described subsequently. Hence, grid resolution was 
not the prime source of the discrepancies. 

The streamwlse velocity profiles across the duct at various downstream 
stations are shown in Fig. 5. As can be seen, there Is fairly good agreement 
between the data and the predictions. However, the axial velocity variation 
along the centerline which is shown in Fig. 6 does not compare as well. In 
Fig. 6 the present calculations are plotted against the experimental data, 
and the predictions of Ref. 3 which uses a different numerical procedure, but 
a similar k-e turbulence model. Both computations indicate a dip in the 
axial velocity which is not observed In the experiments. The precise reasons 
for this behavior is uncertain, but turbulence modeling appears to be the 
most likely culprit. In view of the good agreement with data that was 
obtained for the swirling case, the source of the discrepancies must lie in 
the distinguishing characteristics of the two flows. In reference 4, Brondum 
and Bennet determined the presence of large scale coherent structures in the 
nonswirling case. Furthermore their major effect was precisely in those 
regions which showed the greatest deviation from the experimental data. 
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Therefore, It would appear that the ic-e tranport equations should be 
Investigated to determine how they can be modified to account for the 
phenomena observed In the experiments . 

In Fig. 7 the streamline pattern and In Fig. 8 contours of constant 
axial velocity profiles are shown. The recirculation zone which develops 
along the duct wall extends approximately four duct radii downstream of the 
dump plane, and compares well with the experimental data. 

Since the emphasis of the present effort was to demonstrate the 
capabilities of the numerical scheme, the resolution of the descrepancies 
were not pursued further, but rather the swirling case was considered. As 
initial conditions, the nonswlrl flow field was employed and thereafter the 
swirl velocity was introduced at the Inflow boundary In the annulus over a 
prescribed number of time steps. Two cases were considered at two different 
values of artificial dissipation parameter a. The values of the parameter a 
were .5 and .1, a lower value Indicating less dissipation. 

The results of the computations are shown in Figs. 9 to 14. Both a = .5 
and .1 computations are shown. In Fig. 9 the streamline patterns are shown. 
As can be seen, there is a significant variation In the flow patterns In 
particular In the shape and extent of the recirculation zones. Similar 
differences In the two calcultions can be seen In Figs. 10 and 11 where the 
contours of constant streamwise and azimuthal (swirl) velocities are shown. 

A more dramatic effect is seen in Figs. 12 to 13 where the axial and 
azimuthal velocity profiles are shown at different streamwise locations in 
the duct. The a = .1 calculations are in excellent agreement with the data. 
Compared with the a = .5 solution, the effect of artificial dissipation is 
clearly evident, in the smearing of the profiles and the cutting off of the 
peaks. The final plot figure 14 shows the streamwise velocity distribution 
along the duct. Here again the effect of reduced artificial dissipation is 
evident. In that the computed results are in better agreement with the 
experimental data. 

In comparing the streamwise velocity distributions along the duct. It is 
noted that the velocity at the dump plane is lower than the experimental 
value. Hence, there appears to be a discrepancy with regard to the mass flow 
through the system. Inasmuch as the parameter describing the flow is the 
Reynolds number, the discrepancy in mass flow would vary the Reynolds number 
under which the calculation was run. The reasons for the discrepancies in 
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mass flux are uncertain, but errors In the speclficaton of inflow profiles 
could have played a role. Furthermore, the choice to set function conditions 
at the upstream inflow boundary may have also contributed to the observed 
discrepancies. For subsonic inflow, the preferred inflow boundary condition 
is to set the stagnation pressure, and let the streamwise velocity profile 
adjust to accomodate the mass flux. Since the stagnation pressure was 
unavailable, the mass flux was specified instead. It is felt that further 
investigation is warranted with regard to this aspect of the calcuation. 
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CONCLUSIONS 


Navier~Stok.es calculations employing a k-e turbulence model were 
obtained for the flow resulting from confined swirling and nonswirling 
confined coaxial jets. The calculations compared well with the experimental 
data of Johnson and Bennet and Roback and Johnson, in particular for the 
swirling case. The calculations indicate that, as long as the turbulence 
model is well specified, the numerical procedure can give results that 
compare very well with the data. Furthermore the results demonstrate that 
artificial or numerical dissipation can have a significant effect on the 
computations and must be carefully controlled to obtain accurate predictions. 
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SKETCHES OF TEST SECTION INLET REGION WITH VELOCITY AND COORDINATE SYSTEM 

30 DEG MEAN ANGLE SWIRLER )N ANNULAR INLET DUCT 



Figure 2. Sketch of Test Section Inlet Region (from Ref. 2) 








Figure 3. Schematic of Computational Domain. 
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Figure 4. Computational Grid (Downstream Section not Shown). 
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Figure 5. Axial Velocity Profiles for Nonswirling Flow. 
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Figure 5(a). Axial Velocity Profiles for Nonswirling Flow. 
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Figure 5(b). Axial Velocity Profiles for Nonswirling Flow. 




= 305 mm (exp.), 312.9 mm (pred.) 
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Figure 5(c). Axial Velocity Profiles for Nonswirling Flow 




Figure 6. Mean Axial Velocity Distribution Along Duct 
Centerline Nonswirling Flow. 
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Figure 7. Streamline Pattern for Nonswirling Flow. 
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Figure 12(a). Mean Axial Velocity Profiles for Swirling Flow. 
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Figure 12(b) 
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Figure 13. Mean Azimuthal Velocity Profiles for Swirling Flow. 
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Figure 13(a). Mean Azimuthal Velocity Profiles for Swirling Flow. 




Figure 13(b). Mean Azimuthal Velocity Profiles for Swirling Flow. 
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APPENDIX A 


The governing conservation equations in cylindrical-polar coordinates 
are transformed using the Jacobian transformation, 


y j = y^x,, x 2 , x 3 ,t) 


(A-l) 


T = t 

where (x^ , X 2 , X 3 ) = (r, 0, z). The resulting equations may be written in 
the following compact form: 


a( jw ) 

dr 




i=. L ^ 


(A-2) 


+ 'Jvl * Jyl, i + avT (Jyl 'i G|) } + JS + jc 


dyJ 


dyj 


where 


y,« s 


y J ,j = 


iz! 

at 

dy J 

ax : 


(A-3) 


Further, the coefficients 3£, Y£, ?£ are given by 


0. = -r • 

r. = 1 ■ 

t = JL 

= 1 r™ * 


f • *3 = 


r 2 

C 2 


^3 = 1 


, ' 


(A-4) 


and m = 1 for all equations except the X2~direction momentum equations for 
which m = 2. 
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The vector variables used in Eq. (A-2) are defined as 
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(A-5) 

where n = 1 for i = 1 and n = 0 for i = 2, 3. 


(A- 6 ) 
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(A- 7 ) (A- 8 ) (A-9) 

Note that the velocity componetns (Uj , U 2 , U 3 ) are the cylindrical-polar 
velocity components wntten in cylindrical-polar coordinates. The molecular 
and turbulent stress tensors may be written as 


T ij - T^.. (v - u) Vt< k b v - u -/> k ' S ij 


(A- 10 ) 


and the rate of strain tensor components in cylindrical-polar coordinates are 
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D„ = 


aui 
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- _ i au 2 u, 
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d 33 = 
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(A-ll) 


and 


_ _ i a , . i au 2 . au 3 

v - u = — — — (rUi)+ — + — — 

r 3x, dx 2 ax 3 


(A-12) 


The ' derivatives required in Eqs. (A-ll) and (A-12) must be expressed in terms 

of the computational coordinates yJ using the chain rule. 

• ... 

Finally, the vecotr S contains source terms and certain differential 

terms w£ich do not conform to the basic structure of Eq. (A-2), and the 

vector C contains the additional curvature terras due to the cylindrical-polar 

coordinate system. 
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APPENDIX B - SOLUTION PROCEDURE 
Background 

The solution procedure employs a consistently-split linearized block 
implicit (LBI) algorithm which has been discussed in detail in [5, 6] . 

There are two important elements of this method: 

(1) the use of a noniterative formal time linearization to 
produce a fully-coupled linear multidimensional scheme which 
is written in "block implicit" form; and 

(2) solution of this linearized coupled scheme using a consistent 
"splitting" (ADI scheme) patterned after the Douglas-Cunn Il7 ] 
treatment of scalar ADI schemes. 

The method is thus referred to as a split linearized block implicit (LBI) 
scheme. The method has several attributes: 

(1) the noniterative linearization is efficient; 

(2) the fully-coupled linearized algorithm eliminates Instabilities 
and/or extremely slow convergence rates often attributed to methods 
which employ ad^ hoc decoupling and linearization assumptions to 
identify nonlinear coefficients which are then treated by lag and 
update techniques; 

(3) the splitting or ADI technique produces an efficient algorithm 
which Is stable for large time steps and also provides a means for 
convergence acceleration for further efficiency in computing steady 
solutions; 

(4) intermediate steps of the splitting are consistent with the 
governing equations, and this means that the "physical" boundary 
conditions can be used for the intermediate solutions. Other 
splittings which are inconsistent can have several difficulties In 
satisfying physical boundary conditions [6]. 

(5) the convergence rate and overall efficiency of the algorithm are 
much less sensitive to mesh refinement and redistribution than 
algorithms based on explicit schemes or which employ ad hoc 
decoupling and linearization assumptions. This is important for 
accuracy and for computing turbulent flows with viscous sublayer 
resolution; and 
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(6) the method le general and is specifically designed for the 

complex systems of equations which govern multlscale viscous flow 
In complicated geometries. 

This same algorithm was later considered by Beam and Warming [23], but the 
ADI splitting was derived by approximate factorization Instead of the 
Dougla6-Gunn procedure. They refer to the algorithm as a “delta form” 
approximate factorization scheme. This scheme replaced an earlier non-xielta 
form scheme [24], which has inconsistent intermediate steps. 

Split LBI Algorithm 
Linearization and Time Differencing 

The system of governing equations to be solved consists of three/four 
equations: continuity and two/three components of momentum equation in 

three/four dependent variables: p, u, v, w. Using notation similar to that 

in [5], at a single grid point this system of equations can be written in 
the following form: 

3H($)/3t = D($) + SU) (1) 

where 4> is the column-vector of dependent variables, H and S are column- 
vector algebraic functions of and -D is a column vector whose elements are 
the spatial differential operators which generate all spatial derivatives 
appearing in the governing equation associated with that element. 

The solution procedure is based on the following two-level implicit 
time-difference approximations of (1): 

(H n+1 - H n )/At - B(D n+1 + S n+1 ) ( 1-0) (D n + S n ) (2) 

where, for example, H n "^l denotes H(4> n+ ^) and At = t n+ ^ — t n . The 

parameter 0 (0.5-0- 1) permits a variable time-centering of the scheme, 
with a truncation error of order [At^, (0 - 1/2) At]. 

A local time linearization (Taylor expansion about $ n ) of requisite 
formal accuracy is introduced, and this serves to define a linear differen- 
tial operator L (cf. [5]) such that 

n+1 n n, n+1 n. . . 2, , . . 

D = D + L ( 4> — 4>) + 0( At ) (3) 

Similarly , 

H n+1 = H n + (3H/3<J>) n U n+1 - <|> n ) + 0 (At 2 ) (4) 

S n+1 = S D + (3S/34>) n (4> n+1 -<]>")+ 0 (At 2 ) (5) 
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Eqs. (3-5) are inserted into Eq. (2) to obtain the following system which i6 
linear in <J> n+ l 

(A - 0At L ) (4> - $ ) = At (D + S ) (6) 

and which is termed a linearized block implilcit- (LB1) scheme. Here, A 
denotes a matrix defined by 

A = (3H/3<j>) n - BAt (3S/3<{.) n (7) 

Eq. (6) has 0 (At) accuracy unless H = $* in which case the accuracy is the 
same as Eq . ( 2) . 

Special Treatment of Diffusive Terms 

The time differencing of diffusive terms is modified to accomodate 
cross-derivative terms and also turbulent viscosity and artificial dissipa- 
tion coefficients which depend on the solution variables. Although formal 
linearization of the convection and pressure gradient terms and the resulting 
implicit coupling of variables is critical to the stability and rapid con- 
vergence of the algorithm, this does not appear to be important for the 
turbulent viscosity and artificial dissipation coefficients. Since the 
relationship between p e and dj and the mean flow variables is not conven- 
iently linearized, these diffusive coefficients are evaluated explicitly at 
t n during each time step. Notationally , this .is equivalent to neglecting 
terms proportional to 3 p e / 3$ or 3dj/3$ in L n , which are formally pre- 
sent in the Taylor expansion (2), but retaining all terms proportional to 
p e or dj in both L n and D n . 

It has been found through extensive experience that this has little if 
any effect on the performance of the algorithm. This treatment also has the 
added benefit that the turbulence model equations can be decoupled from the 
system of mean flow equations by an appropriate matrix partitioning (cf. 

[6]) and solved separately in each step of the ADI solution procedure. This 
reduces the block size of the block tridiagonal systems which must be solved 
in each step and thus reduces the computational labor.' 

In addition, the viscous terms in the present formulation include a 
number of spatial cross-derivative terms. Although it Is possible to treat 
cross-derivative terms implicitly within the ADI treatment which follows, it 
is not at all convenient to do so; and consequently, all cross— derivative 
terms are evaluated explicitly at t n . For a scalar model equation 
representing combined convection and diffusion. It ha6 been shown by Beam and 
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Warming 1 2 5 J that the explicit treatment of cross-derivative terms does not 
degrade the unconditional stability of the present algorithm. To preserve 
notations! simplicity, it is understood that all cross— derivative terms 
appearing in L n are neglected but are retained in D 11 . It is Important to 
note, that neglecting terms in L n has no effect on steady solutions of Eq. (7), 
since - 4> n = 0, and thu6 Eq. (7) reduces to the 6teady form of the 

equations: D n + S n ■* 0. Aside from stability considerations, the only 

effort of neglecting terms in L n is to introduce an 0 (At) truncation error. 

Consistent Splitting of the LBI Scheme 

To obtain an efficient algorithm, the linearized system (7) i6 split using 
ADI techniques. To obtain the split scheme, the multidimensional operator L is 
rewritten as the sum of three "one-dimensional “ sub-operators (i «= 1, 2, 3) 
each of which contains all terms having derivatives with respect to the i-th 
coordinate. The split form of Eq. (7) can be derived either as in [5, 6] by 
following the procedure described by Douglas and Gunn [17] in their 
generalization and unification of scalar ADI schemes, or using approximate 
factorization. For the present system of equations, the split algorithm is given 
by 

(A - 0AtL") (t* - <j. n ) = At <D n + S n ) ( 8a ) 

n A n 

(A - BAtL") («}> - O = A (* - O ( 8b) 

(A - BAtL") (* n+1 - <p n ) = A <♦** - <p U ) ( 8c) 

where $>* and $** are consistent intermediate solutions. If spatial 
derivatives appearing in and D are replace by three-point difference 
formulas, as indicated previously, then each step in Eqs. (lOa-c) can be solved 
by a block-tridiagonal elimination. 

Combining Eqs. ( 8 a -c) gives 

(A - BAtL") A _1 (A - BAtL^) A _1 (A - BAtL") (<f> n+1 - <p n ) = At (D° + S°) ( 9 ) 

which approximates the unsplit scheme (8) to 0 (At^). Since the intermediate 

steps are also consistent approximations for Eq. (8), physical boundary 
conditions can be used for <f>* and ■}>** [5, 6]. Finally, since the 

are homogeneous operators, it follows from Eqs. ( 8a— c) that steady solutions 
have the property that = <j>* = <|>** = and satisfy 

D n + S n = 0 (10 ) 

The steady solution thus depends only on the spatial difference approximations 
used for (10), and does not depend on the solution algorithm itself. 
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16 Abstract 


The existence of large-scale coherent structures In turbulent shear flows has 
been well documented. Discrepancies between experimental and computational data 
suggest a necessity to understand the roles they play In mass and momentum trans- 
port. Using conditional sampling and averaging on coincident two-component 
velocity and concentration-velocity experimental data for swirling and nonswirl- 


ing coaxial jets, triggers for Identifying the structures were examined. Concen- 
tration fluctuation was found to be an adequate trigger or Indicator for the 
concentration-velocity data, but no suitable detector was located for the two- 
component velocity data. The large-scale structures are found In the region 
where the largest discrepancies exist between model and experiment. The tradi- 
tional gradient transport model does not fit In this region as a result of these 
structures. The large-scale motion was found to be responsible for a large per- 
centage of the axial mass transport. The large-scale structures were found to 
convect downstream at approximately the mean velocity of the overall flow In the 
axial direction. The radial mean velocity of the structures was found to be 
substantially greater than that of the overall flow. 
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